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Abstract. Wc set up a methodology for computer assisted proofs of the existence and the KAM 
stability of an arbitrary periodic orbit for Hamiltonian systems. We give two examples of application for 
systems with 2 and 3 degrees of freedom. The first example verifies the existence of tiny elliptic islands 
inside large chaotic domains for a quartic potential. In the 3-body problem we prove the KAM stability 
of the well-known figure eight orbit and two selected orbits of the so called family of rotating Eights. 
Some additional theoretical and numerical information is also given for the dynamics of both examples. 



1. Introduction 

KAM Theorem (see [14l[2][2T| and also [3]) is a fundamental result for Hamiltonian systems because it 
ensures the existence of a set, nowhere dense but of positive measure, of points of the phase space which 
behave in a regular, quasi-periodic way. The main point is that the system should be a perturbation of 
an integrable system and that a non-degeneracy condition, asking for the invertibility of the actions to 
frequencies map, has to be satisfied. The standard notation, being a Hamiltonian system given by 

(1) q = H p p = -H q , 

where the Hamiltonian H(q,p) : — » R is a smooth function defined on an open set 57 C R™ +1 x R™ +1 , 
will be used. 

If we consider the dynamics close to a fixed point the methodology is simple. Assume that the fixed 
point is totally elliptic or the problem can be reduced to the totally elliptic case, for instance by restricting 
the attention to the centre manifold. Then one can proceed to compute the normal form up to a moderate 
order, say to order 4 in the (q,p) variables. Assuming that no resonances occur up to this order then one 
can consider the normal form as the integrable Hamiltonian and the remainder as the perturbation, and 
it is easy to check the non-degeneracy condition. This approach has been used, e.g., in the study of the 
vicinity of the collinear libration points in the general planar three-body problem, restricted to the centre 
manifold, see [17] . A moderate number of arithmetic operations allows to decide on the applicability of 
KAM Theorem. 

The problem is much more involved when we want to apply KAM Theorem around an arbitrary totally 
elliptic periodic orbit which is not known analytically. Even if some analytic expression of the orbit is 
available, the study of the dynamics on the vicinity at the required order can be not feasible analytically. 
As it is usual, one can restrict the problem to the study of the vicinity of a fixed point of a symplectic 
map on a suitable Poincare section in dimension 2n. 

The goal of this paper is to set up a methodology for the rigorous check of the KAM conditions for 
the symplectic map (see, e.g., [3]). 

We give two examples of application. The first one is a simple classical Hamiltonian system with two 
degrees of freedom and depending on a parameter c. The main feature is that the potential consists only 
on quartic terms. Changing c the system can be integrable or display large chaotic domains. In these 
domains one can guess, by numerical computation of iterates of a Poincare map, that some tiny islands 
exist. The problem is to show, rigorously, that indeed there are elliptic periodic points of the Poincare 
map inside these islands and KAM conditions hold. 

The second example concerns the well-known figure eight solution of the general three-body problem 
with equal masses. See [7] for a proof of the existence of that orbit, found numerically by Moore [T5]. This 
is an example of "choreography" (see |25j . where the notion of choreographic solution was introduced, 
and the references therein), that is, a T-periodic solution of the TV-body problem (TV = 3 in the present 
case) such that all the bodies move along the same path with time shift T/N between consecutive bodies. 
This topic for N = 3 has been studied by present authors in [10] where, in particular, it was proved 
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the totally elliptic character of the figure eight on fixed energy levels and remaining at the zero level of 
angular momentum. Using reductions the problem becomes a Hamiltonian with three degrees of freedom. 
The related Poincare map is 4D. In |24j it was claimed, based on a non-rigorous high order computation 
of a normal form, that the KAM condition is satisfied around the figure eight. The computation of the 
local expansion of the Poincare map was done by numerical differentiation using multiple precision and 
optimal step size for the different orders. In the present paper the validity of the KAM condition for the 
figure eight orbit is established rigorously. 

The paper is organised as follows. In the next two sections the examples with 2 and 3 degrees of 
freedom are presented and several relevant properties of them are proved or mentioned. In particular the 
reduction of the three-body problem in present case is explicitly carried out, based on |27| . Then the 
methodology to be applied is explained, introducing the required notation and emphasizing the rigorous 
aspects of the CAP (Computer Assisted Proofs). Finally the results obtained by applying the methods 
to both examples are shown. 



H= -(p 2 x +pl + x 4 + cx 2 y 2 +y 4 ) 



2. A FAMILY OF QUARTIC POTENTIALS 

As a first example we consider a very simple Hamiltonian 

1 

2 ( 

where c > —2 is a real parameter. This is a system widely considered as a paradigm of chaotic system for 
c large in the relations between classical and quantum mechanics, see for instance [3J [5] and references 
therein. The energy should be positive and, due to the homogeneity, it can be considered equal to a fixed 
value. We shall consider the level H = i. Note that for c < —2 unbounded motion occurs. The system 
has some obvious symmetries: It is reversible with respect to time and the changes of sign of x and/or 
y leave the equations invariant. Furthermore, the symplectic change induced by the change of variables 
(x, y) (u,v) = -^(x + y, y — x) keeps the form of the Hamiltonian with the parameter c replaced by 
12 - 2c 

c = , after a scaling to normalise the coefficients to x and y to 1. Obviously the map c — > c 

2 + c 

16 

is an involution having c = 2 as fixed point. It can be written also as c + 2 = — - . When c ranges in 

(—2, oo ) increasing its value, the parameter c ranges in the same interval but decreasing. 

It is immediate to check that the planes (x,p x ), (y,p y ), (x = y,p x = P y ) and (x = —y, p x = —p y ) are 
invariant and, for the last two, modulo the change c f> c, we have the same phase portrait than for the 
first two. 

The first question to be addressed is the integrability of the Hamiltonian. To this end we observe that 
H = T + V is a classical Hamiltonian and V is homogeneous of degree k = 4. We specialize to this degree 
of homogeneity a theorem due to Morales-Ramis [20] 

Theorem 1. Assume H = T + V = h(p1 + ... +P„) + V(xi, . . . , x n ) where V is homogeneous of degree 
4-. Let z — (zi, . . . , z n ) be a solution of z — \7V(z) and let Ai, . . . , A„ be the eigenvalues of Hess V(z). 
Then, if H is completely integrable with meromorphic first integrals, the values of A; must be equal to 
numbers of the form 2m 2 — m, 2m 2 + |m + X or 2m 2 + 2m + | for m E Z. 

Corollary 1. The Hamiltonian H is only integrable for c = 0, 2,6. 

Proof. Beyond the trivial eigenvalue Ai = 3, due to the homogeneity, one has A2 = c/2, which should 
be of one of the forms above. Also c/2 should be of one of these forms. This reduces the possible values 
of c to 0, 2, 6. On the other hand the case c = is obviously separable and it is also c = which 
corresponds to c = 6. Finally in the case c = c = 2 the symplectic change x = r sirup, y = rcosip, p x = 
p v cos ip/r + p r sin ip, p y = —p v sin ip/r + p r cos ip converts the Hamiltonian to H = ^ (p 2 /r 2 + p 2 , + ?* 4 ) , 
reducible to 1 degree of freedom. □ 

See Figure Q] for an illustration of the phase portrait of V for values of c close to the integrable cases. 
In the left (resp. right) column the value of c is smaller (resp. larger) than the integrable one. The plots 
allow to identify easily the bifurcations which occur at the integrable cases. 

To understand the dynamics of H as a function of c it is useful to consider the Poincare section S 
through y = 0,p y > 0, defined for p 2 + x A < 1. The boundary of that domain is a periodic orbit on 
the invariant plane y = p y = 0. The initial data x = p x = give rise to a periodic orbit y = y(t), 
sitting on the (y,p y ) plane, which corresponds to a fixed point of the Poincare map V . The normal 
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Figure 1. Phase portrait on S, the Poincare section y = 0, close to the integrable cases. 
Left column from top to bottom: c = —0.1, 1.9, 5.9. Right column c = 0.1, 2.1, 5.9. The 
variables displayed are (x,p x ). 



variational equations £' = 77, rf = — cy(i) 2 £ give the linear stability of (0,0). It is easy to check that 
the trace of DV at the origin decreases from +00 to 2 for c < and then, for increasing c, it oscillates 
between -2 and 6. In particular it takes the value +2 for c = Ck — k(k + l),fc e N. For the values 
c = (k — l/2)(fc + 1/2) it takes, alternatively, the values +6 and —2 for k even and odd, respectively. 
See Figure [2] for an illustration of the behavior of the trace. The set of values of c for which |Tr| < 2 in 
the union of the intervals (co, Ci) U (02, C3) U . . ., where one has linear stability. The first linear stability 
intervals are (0,2), (6, 12), (20,30), .... It is also clear that the periodic orbit at the boundary of E has 
the same stability properties as the orbit we have just discussed. 

In a similar way one can consider the initial conditions x = 0, p x = ^/l/2 which correspond to a 
periodic orbit in x = y,p x — Py In that case one can take x + y = as Poincare section. Replacing 
c by c one has a similar result: For c — > — 2 + , which implies c — > 00, this fixed point is unstable for 
c > 6. Stability intervals for c are identical to the ones given before for c in the case of the fixed 
point at the origin. They give rise to intervals which accumulate to —2. The first intervals (in the c 
parameter) are (2,6), (-6/7,0), (-3/2,-14/11), (-50/29,-18/11), (-42/23), (-66/37), etc. We shall 
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Figure 2. Trace of the differential of the Poincare map V at the origin as a function of 
c. On the vertical axis, to prevent from the effect of the large values when c — > — 2 + , we 
plot argsinh(Tr) instead of Tr. The horizontal lines correspond to —2, 2 and 6. 

refer to these two periodic orbits as the basic ones. Of course, symmetries give rise to similar orbits, like 
the one through x = 0, p x = —\J\j2 

To see the evolution of the phase space away from the integrable cases we have computed an estimate 
of the "fraction of chaotic motion" in £ as a function of c. Due to the symmetries it is enough to do 
the computations for x > 0, p x > 0. In the domain bounded by x = 0, p x = and p x + x 4 = 1 we have 
selected "pixels" with centre of the form (i/2000, j/2000), i,jgN and for each one we have estimated the 
maximal Lyapunov exponent A. In fact we are not interested on the concrete value but rather on whether 
one can accept A = 0. Symmetry and some other simplifications allow to reduce the computational task. 
Then the fraction of chaotic motion tp(c) is estimated as the number of pixels for which one has evidence 
that A > divided by the total number of pixels in the domain. Several checks have been done using 
different strategies and maximal number of iterates of the Poincare map to have reliable information (see, 
e.g., [23] for details). The results are shown in Figure [3] 
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Figure 3. Fraction of chaotic dynamics 4>(c) m the Poincare section £ as a function of 
the parameter c. To better see the behavior of the function tp the variable displayed in 
the horizontal axis is d(c) = log(c + 2)/log(2) instead of c. The three integrable cases 
c = 0, 2, 6 for which ip(c) = are seen on d = 1, 2 and 3, respectively. 

The interpretation of the figure is clear: close to the integrable cases most of the dynamics is regular. 
In fact, if c* is one of the values of c for which one has integrability, the fraction of chaoticity seems to 
be exponentially small in |c — c*| for nearby values. Between these values of c the value of tp( c ) is below 
0.12. Then, to the left of c = and to the right of c = 6 there is a quick increase of tp(c). But at the 
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ranges in which one of the basic periodic orbits is linearly stable, we can expect the existence of islands, 
which decrease the measure of the chaotic domain. The oscillations of the decrease are becoming smaller 
when the limits, either c — > — 2 + or c — > oo are approached. 

More concretely, the system away from the range of values of c where it is close to integrable, has a 
well defined and repetitive structure. To this end we consider the "fraction of integrability" 1 — ip( c )- 
Figure |4] shows it for the values of the parameter c on the right hand side of the domain which contains 
the close to integrable dynamics. On the left hand side of that domain the plot is symmetrical to the one 
shown in Figure 0] left. To better sec the scaling properties, the horizontal parameter r is an extension 
of the index k to the real numbers, defined by r(r + 1) = c. The different "bumps" are quite similar, 
and the heights scale as « 1/r. The right part of Figure H] shows the behavior for c 6 [C4, ce], which is 
similar for all the ranges [c2 m , C2m+2\, m > 2. At c = 20 the central periodic orbit becomes again elliptic 
(see Figure [2]) (remember that the same thing happens for the periodic orbit at the boundary of £). A 
domain of regular motion starts which is reminiscent of the behavior found in many other problems. See, 
for instance, the reference |26) in the case of the Henon map. However there are important differences 
due to the many symmetries of present problem. Anyway the mechanisms to explain the behavior of the 
plot are the ones explained in [26]. For instance, around c = 25.6 the "last" curve surrounding the period 
4 islands breaks down and this increases the size of the main chaotic zone. 

At c = 30 the central periodic orbit becomes unstable by a pitchfork bifurcation and its separatrices 
display the typical figure eight shape, similar to the example displayed in the central plot in Figure [5J 
Each part of the figure eight, skipping the invariant curves surrounding the full figure eight scparatrix, is 
really close to the behavior of the Henon map (except by nonlinear scalings in variables and parameters). 
The decrease of the regular dynamics in Figure 0] right from c = 30 to c « 31.1 is due to the destruction 
of the invariant curves surrounding the full figure eight scparatrix. Then, at the approximate c values 
32.0, 32.5, 33.5, 35.0 and 36.8 the decrease in the fraction is due to the destruction of curves surrounding 
period 6, 5, 4 and 3 islands and to the period doubling, respectively. It is worth to mention that for the 
case of period 3 almost no point has regular dynamics. The dynamics in that case, including existence 
of tiny period 3 islands and other properties, is in perfect agreement with the ones of the Henon map as 
studied in [TTj . 




4 8 12 16 20 25 30 35 40 

Figure 4. Left: The fraction of integrability as a function of r, defined by r(r + 1) = c. 
Right: Idem for a range of c containing [C4, ce] = [20,42] using c as horizontal variable. 



Our goal is to detect periodic orbits, well inside a chaotic domain, to which we can apply the method- 
ology to be explained to test the KAM conditions. To this end we have selected the value c = —0.90 (or, 
equivalently c = 138/11), for which the value of ip is approximately 0.9528. Iterates of the Poincare map 
are easily obtained. Using Taylor method at order 30 and a maximal relative truncation error of 10 -21 
they are computed at an average rate of 15, 000 iterates per second for that value of c. 

The Figure [5] shows a global view of the phase portrait and some details. As this value of c is in one 
of the instability domains, the periodic orbit through (0, \/l/2) is unstable, but c is not too far from 
the left end of one of the stability intervals, at c = —6/7. Then the chaotic zone around the figure 
eight separatrix is still surrounded by invariant curves. Inside that chaotic zone there exist several stable 
periodic orbits. On the central part of Figure [S] one can see a period 16 orbit. We shall test the KAM 
conditions for that orbit. An approximate initial condition on X is x = 0, p x = 0.51780665799545. But 
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Figure 5. Phase portrait on £ for c = —0.90. Left: a global view (in grey) and two 
periodic orbits (in black). Centre: a magnification of the chaotic zone around the figure 
eight separatrix associated to the unstable periodic orbit on (0, y/T/2) and a stable 
periodic orbit of period 16. Right: The island and some satellite islands near one of the 
period 3 points, well inside the large chaotic domain. 



also a periodic orbit of period 3, far away from that separatrix, deeply inside the chaotic domain, has been 
found. An approximate initial condition for this orbits is x = 0.352718557335, p x = 0.547882838499, and 
the applicability of KAM theorem for this orbit will also be checked. 

We remark that the largest area inside an invariant curve around that pcriod-3 island is of the order 
of 7 x 10~ 6 . It is not easy to capture this periodic orbit, but the previous computation of Lyapunov 
exponents in a grid with small stepsize is of great help. 



3. The problem of the KAM stability of the Eight 

The figure eight periodic orbit, shown in Figure [6j is a remarkable solution of the planar three-body 
problem with equal masses [7]. The three bodies move on the plane along the same path in solutions of 
the form q(t), q(t + T/3), q(t + 2T/3) where T is the period. In [24] a detailed numerical study of this 
orbit and an extended vicinity of it was done, looking also at the effect of small changes in the masses 
and the bifurcations that they create. See also for choreographies related to the figure eight, like the 
satellite and the relative ones. The numerical evidence given in |24j suggested that non only the orbit 
was linearly stable but also KAM theorem applies around it. The rigorous proof of the linear stability 
was given in [TU] and the proof of the applicability of KAM is studied now. 




The figure eight solution, three nearby partially hyperbolic periodic orbits, as well as several 2D and 
3D tori around the eight can be seen in |24j . The orbit has zero total angular momentum and due to the 
homogeneity of the potential and Kepler's law, one can fix either the level of (negative) energy or the 
period. 

The first step to be done is the reduction of the problem to a three degrees of freedom system. This 
is a classical result and we follow the exposition that can be found in Whittaker's treatise [27]. For 
completeness a short account is given below. 
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3.1. Reduction of the 3-Body Problem. We shall assume that the masses are equal to 1. Let 
Pa j Cl2, q 3 be the positions of the three bodies in R 2 and Pi,P2,P 3 the corresponding momenta. The 
Hamiltonian is 

H = i(| P i| 2 + | P2 | 2 + |p 3 | 2 ) - . 1 ; - ; 1 . - . 1 | 

and the angular momentum is M = qi A pi + q2 A P2 + q 3 A p 3 . 

Let qi, q 2 , q 3 , p[, p 2 , and p' 3 be new variables introduced by means of the generating function 

W = (Pi,qi) + (P2,q 2 ) + (Pi + P2 +p 3 ,q 3 ), 
where ( , ) denotes scalar product. We recall that 

dW . dW 
q ^9p7' ^ = W- 3 = W 

and, hence, the change gives 

qi = q' : + q 3 , q 2 = q' 2 + q^ q 3 = q 3 , Pi = pi, P2 = P 2 , P 3 = P 3 - Pi P 2 - 

Because of the centre of mass integrals, it is not restrictive to assume qi + q 2 + q 3 = 0, pi + p 2 + p 3 = 0, 
which amounts to q 3 = — (q^ + q' 2 )/3 and p' 3 = 0, Hence, the new expressions of H and M, skipping 
the ' for simplicity, are 

H = |pi| 2 + |p 2 | 2 + (Pi,p 2 ) - Ar - 7— T - 1 " 1' M = Qi Api +q 2 Ap 2 , 

|qil [qui |qi-Q2l 

which reduce the system to 4 degrees of freedom. 

Let (91,(72) be the components of qi, ((73,(74) the ones of q 2 and, in a similar way, we define the 
components of the p variables. We introduce the generating function 

W = piq[ cos q' 4 + p 2 q[ sin q 4 + p 3 (q' 2 cos q' 4 - q' 3 sin q' 4 ) + Pi{q' 2 sin q' 4 + q' 3 cos q 4 ), 
which gives raise to the transformation 

qi = q[ cos q' A , q 2 = q[ sin q' A , q 3 = q' 2 cos q 4 - q' 3 sin q 4 , q 4 = q' 2 sin q' 4 + q' 3 cos q 4 , 
pi = p[ cos 94 - ((P4 - q 2 p 3 + q 3 p 2 ) sin q 4 )/q[, 
p 2 = p\ sinq 4 + ((P4 - ^3 + q 3 p 2 ) cos q 4 )/q[, 
P3 = P 2 cos q' 4 - p 3 sin q 4 , p 4 = p' 2 sin q 4 + p 3 cos q 4 . 
Skipping again the ', one can write the Hamiltonian and angular momentum in the new variables. It 
turns out that the new q^ does not appear in H. It is a cyclic variable. Hence, the conjugated variable pi 
is constant. But it is immediate to see that M = pa- We still keep its value in the Hamiltonian, despite 
in the case of the figure eight M = 0, because it plays a role in the "rotating eights" solutions. The final 
reduced Hamiltonian has the form 

(2) H = p\ +pl +pj +pip 2 - — {q 2 p 3 - q-3P2 - Pa) + \{q2Pi - qzPi - Pa) 2 

qi <7i 

1 1 1 



91 (<Z 2 2 + <7 3 2 ) 1/2 ((?i-92) 2 + ? 3 2 ) 1/2 ' 
The new variables have a simple geometrical meaning. Let us denote the positions of the masses 
as m,j,j = 1,2,3 and as mirrij the vector from m, to mj. Then q 4 is the norm of m 3 m\ and p\ is 
the component of the linear momentum of m\ projected along m 3 mi; q 2 and q 3 are the projections of 
m 3 m 2 along m^mi and orthogonal to it and, in a similar way, p 2 and p 3 are the projections of the linear 
momentum of m 2 ; finally (74 is the angle between the x-axis and 7713ml. As said, pa = M. 

3.2. Rotating Eights. If the angular momentum M goes away from zero, the periodic orbit can be 
continued. It becomes quasi-periodic with two basic frequencies, which will produce a periodic orbit if 
the ratio of frequencies is rational. But it can be seen again as a periodic solution, even as a choreography, 
using a rotating frame, a fact noticed by M. Henon [9j who also found that the continuation leads to 
collision orbits. Due to the symmetries of the problem it is enough to consider M increasing to positive 
values. For most of these "rotating eights" we can apply the same algorithm as for the Eight. The Figure 
[7] shows the values of ol\ < a 2 such that the eigenvalues of the Poincare map associated to the rotating 
eight are exp(±27riai), exp(±27ria2), as a function of M. The monotonically decreasing line gives the 
minimal distance between the bodies along the orbit, an evidence of a nearby collision. The plot changes 
at M = M p d w 0.4467, when a 2 = 1/2, and then the orbit loses stability with the appearition of a period 
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doubling bifurcation. For M > M p d we plot the logarithm of the modulus of the dominant eigenvalue 
(divided by 2 to fit in the plot). Low order resonances k\oti + k^ct2 £ N with |fci| + \k?\ < 4 are easy 
to detect for different values of M. For them additional terms will appear in the normal form making 
our arguments not valid in these cases. Therefore, as an example, we have selected two values of the 
angular momentum, Mi = 0.0048828125 and Mi = 0.1484375, that are far enough from resonances. A 
non-rigorous exploration of orbits close to the rotating eights for \M\ < M p d gives evidence that for most 
of them the numerical simulation suggests the existence of tori, but for some resonances they seem to be 
destroyed. 




0.1 0.2 0.3 0.4 0.5 0.6 

Figure 7. Arguments a±, of the eigenvalues cxp(±27ri exp(±27ri a^) correspond- 
ing to the rotating Eights solutions for a range of values of angular momentum M. When 
«2 reaches the value 1/2 for M = M p d ~ 0.4467 a period doubling is produced and the 
orbit becomes unstable with a dominant eigenvalue A < — 1. From that value on, instead 
of «2 we plot log(|A|)/2. The line which decreases monotonically shows the minimal 
distance between the bodies along the rotating Eight. 

Initially the rotating eights, in coordinates which rotate with the suitable angular velocity, look similar 
to the orbit shown if Figure |H1 but the left and right hand side lobes are no longer symmetric. Later 
on the orbits in the rotating frame can develop extra loops. As an example Figure |S] shows the solution 
obtained for M k 0.4493, shortly after the orbits become unstable. The value of M has been selected 
so that in the fixed frame (middle panel) the orbit is also periodic. On the left panel the orbit is shown 
in a rotating frame. At some moment, as displayed, one of the bodies is on the rightmost point on the 
path, on the small loop on the right, and the other two are symmetrical w.r.t. the horizontal axis, on 
the large loop on the left. When the bodies move they pass very close to collision on the tiny central 
lobe. The right panel shows a 3D projection, on the variables (q±, (72, Pi), of what is seen in the Poincarc 
section. The large dot represents the periodic orbit. An initial point, taken by adding 10~ 9 to qx, is 
iterated under the Poincare map. The points are scattered close to the "separatrix" associated to the 
unstable/stable directions. Of course, these ID manifolds do not coincide, but the splitting must be 
expected to be exponentially small in M — M p d- 
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Figure 8. A rotating eight for M ~ 0.4493, already on the linear instability domain. 
See the text for a detailed explanation. 
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4. Preliminaries and the algorithm 



Before going into details we want to emphasize that the algorithm presented in this paper is rigorous. 
By rigorous we mean that during all the computations we take into account and bound all possible errors. 
In this way we get not the exact values but verified estimates of the computed quantities. Therefore the 
theorems that we apply have assumptions of a special kind (i.e. inequalities, inclusions etc.) that can be 
checked using those estimates. As a result we obtain a computer assisted method to have proofs of the 
existence and the KAM stability of periodic orbits for Hamiltonian systems that has full mathematical 
rigor. In theory all calculations could be done by paper and pencil, but in practice the number of 
operations, even if they are very basic and trivial, exceeds human resources. 

4.1. Interval arithmetics. As the precision of the computer is finite we use an interval arithmetic 
to take care of the round-off errors. All floating point operations are replaced by the corresponding 
operations on closed intervals such that we always obtain some representable superset of the true result. 
This is also extended to all elementary functions. 

For the rest of the paper, following [13] , we use boldface to denote intervals and objects with interval 
coefficients. For those objects we use the names interval vector (or 602;), interval matrix etc. to stress 
their interval nature. For a set S by hull S we denote the interval hull of S i.e. the smallest product of 
intervals containing S. For an interval a — [a, a] we define its diameter by diamet := a — a. We define 
the diameter of a box x = (x\, x%, . . . , x n ) as a maximal diameter of its components i.e. diam(a;) = 
max{diamxfc|fc = 1, . . . , n}. 

4.2. The rigorous computation of Taylor expansions of Poincare maps. To rigorously integrate 
ODEs and to obtain verified enclosures for the partial derivatives with respect to the initial conditions we 
use C 1 -Lohner [33] and C r -Lohner [25J algorithms implemented in the CAPD library [5]. Those algorithms 
are based on the Taylor integrator and set representations proposed by Lohner. For some set of data 
(say the box a;) and time step h the C-Lohner algorithm produces rigorous enclosures for the solution 

<p(x, h) of the ODE and its derivatives - for all initial points x 6 x and all multiindices \a\ < r. 

The specialized C^-Lohncr algorithm is able to compute only first order derivatives. 

The "naive" representation of the derivatives as a table of interval vectors leads to a huge overestimation 
due to the "wrapping effect" (see [US [16]). Hence, internally during the integration, for each multiindex 

a the corresponding vector v = - — g^' h ^ is stored in one of the Lohner representations |16j . For that 
reason even C° algorithms need C 1 information to properly set coordinates to suppress the "wrapping" 
error. In the current implementation to store derivatives we use doubletons 



where Xq is a point vector (the centre), C and B are matrices of "good" coordinates (usually C is close 
to the Jacobian matrix and B is its orthogonalization) , r represents the initial size of vector v and r 
stores all computational errors. 



equations. If they would be solved directly it would lead to the integration in a high dimensional space 
and is usually inefficient (most of the rigorous solvers internally need C 1 information that squares the 
dimension) . The C-Lohner algorithm makes use of the special structure of the variational equations to 
avoid this and as a result it can bound derivatives to an arbitrary order r in an efficient way. 

On top of those rigorous ODE solvers the CAPD library implements algorithms to rigorously compute 
Poincare maps and their derivatives with respect to the initial conditions for an affine sections (for details 



4.3. Notation and definitions related to Hamiltonian systems. We denote by J the Poisson matrix 



where /„ denote the n dimensional identity matrix. The Poisson bracket of functions /, g : R 2 ™ D M \-t R 
is a new function {/, g} := (V/) T JVg. 

By N + we denote set of all positive integers i.e. N + = {1,2,...} and we define N := {0} U N + . 
An element of N n will be called a multiindex. For a multiindex k = (fc 1; fc 2 , . . . , k n ) and a vector v = 
(yi, v 2 , ■ ■ ■ , v n ) we define 



v = x + Cr + Br, 



To compute all the derivatives up to order r for an n dimensional 




see [2J[28]). 
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• V k = v kl v% 2 ...V, 



1 1 



,/.■ 



A vector (Ai, A2, . . . , A„) £ C™ satisfies the non-resonant condition up to order r if for all multiindices 
k = {k\, &2, . . . , k n ) such that 1 < \k\ < r and all j £ {1, 2, n} we have 



4.4. The algorithm. For a given Hamiltonian system (JTJ we assume that we have an approximate 
initial condition Xq of a periodic orbit and a Poincare section E that contains xo- For fixed level of energy 
ft. « H(xq) the Poincare map defines a symplectic map 



We will use x = (q,p) G R" x M" to denote local canonically conjugated variables. Let xq correspond to 
xq in these variables. 

In this setting to prove existence and KAM stability of a periodic orbit given by approximate initial 
conditions xq it is enough to prove for some k £ N + the existence and KAM stability of the unique fixed 
point of V k in some small neighbourhood of xq . The algorithm consists of the following steps: 

(1) Proof of the existence of the unique fixed point (a periodic orbit). 
Rigorous estimates of initial conditions. 

(2) Proof of the linear stability. 

(3) Computation of a rigorous Birkhoff normal form. 

(4) Checking an appropriate non-degeneracy condition. 

The details for each step will be given in the following subsections. 

4.5. Proof of the existence. In the first step of the algorithm we prove the existence of an unique 
periodic orbit close to xq and obtain rigorous bounds for its initial conditions. Therefore the preliminary 
step is to reduce the Hamiltonian system ((T|) by suitable symplectic transformations so that for a given 
energy level the periodic orbit is isolated. Section [3] contains a (very classical) example showing how it 
was done in the case of the 3-body problem. 

For a proof we take a box x C K 2 ™ with centre in xq , and compute rigorous estimates of the interval 
Newton operator 



for F(x) = V (x) — x. If we succeed to verify that N(xq,x, F) C x then the interval Newton theorem 
[TJ [55] ensures that inside x there exists a unique fc-periodic point of V . Instead of Newton method 
one can also use the interval Krawczyk method |15[ 112] which do not requires the whole interval matrix 
hu\\(DF(x)) to be invertible. 

Remark 1. The problem of proving the existence of zeros of F when, as in present case, involves the 
computation ofV k (x) is well suited for the use of the parallel shooting method. In our implementation 
we make use of it to improve precision and to speed up computations. 

Remark 2. One can improve rigorous estimates of the initial condition of the periodic point by further 
iteration of the interval Newton or Krawczyk operator. 

4.6. Proof of linear stability. Current step goal is to prove that all the eigenvalues (Ai, A2, . . . , A2n) 
of the differential of the iterated Poincare map 

DP = > wh e r e x is a k periodic point of V, lie on 

the unit circle. We want also to obtain rigorous estimates Xj such that Xj £ Xj for j = 1, . . . , 2n. 

The point x is not known exactly, but from the previous step we have rigorous estimates x of it. From 
estimates DP = — dx using e.g. verified root finding methods applied to the characteristic polynomial 
one can obtain estimates of the eigenvalues. Because those estimates are given by some boxes, part of 
them are out of the unit circle. But still the proof of linear stability is possible due to the fact that our 
system is Hamiltonian. 

Lemma 1. Let A £ l 2nx2n be a symplectic matrix with eigenvalues (Ai, . . . , \2n), cmd let Xj be boxes 
such that Xj £ Xj for j = 1, . . . , In. If the following holds 



(3) 



A, / A fc . 



7V(x , x, F) = xo — hull(D^(a;))- 1 F(io) 



(AI) i Im (A,) for j = 1, . . . , 2n 
(A2) (A^nAfc^) =► j = k, 



then all eigenvalues of A are distinct and lie on the unit circle. 
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Proof: The matrix A is symplectic, hence if A is an eigenvalue of A then also Xj = A -1 and A& = A are 
eigenvalues. But then assumptions (Al) and (A2) ensure that A -1 = A and hence ||A|| = 1. From (A2) 
we have also that all eigenvalues arc distinct. □ 
This general method requires sharp bounds for the eigenvalues. This can be a not so easy task in 
general. Another possibility is to translate constraints to the characteristic polynomial. In this case one 
proves first that the eigenvalues are on the unit circle and then rigorous enclosures Aj are obtained using 
this fact (see for example [Tu]). 

4.7. Computation of Birkhoff normal form. The literature on how to compute Birkhoff normal form 
is very rich. There are also general software packages that can do it in a non rigorous way. Here we want 
to explain how to make this process rigorous. All the computations are done using interval arithmetics 
in the way that at the end the normal form will have interval coefficients that contain the exact values. 
Let T(x) = (EfeeN- 

Ci,k% 7 • ■ • iX/fceN 2 " c 2n,kX ) be the Taylor series of an analytic symplectic map 
around a totally elliptic fixed point and let (Ai, . . . , A„, Ai, . . . , A n ) be the eigenvalues of the linear part 
of T. Then for j 6 {1, 2, 2n} and a multiindcx k such that kj = kj± n + 1 and k m = k m +n for m ^ j 
the condition ((3]) is not satisfied and a resonance occurs. We call it an unavoidable resonance and we say 
that the term Cj t kX k corresponds to that resonance. A Taylor series T is said to be non-resonant if only 
unavoidable resonances are present. 

The goal of this section is to make a symplectic change of variables such that in the new variables the 
Taylor series T up to a given order r reads 

z I— > Az + T r (z), 

where A is a diagonal matrix, A = diag{\\, . . . , X n , Ai, . . . , A„), and T r (z) contains only terms corre- 
sponding to unavoidable resonances. This is the so called non-resonant Birkhoff normal form. In what 
follows we present an algorithm that computes the non-resonant Birkhoff normal form up to order 3, but 
it can be easily extended to any given order. 

First, using C r -Lohner algorithm [28] we compute rigorous enclosures of the coefficients of the Taylor 
expansion of V k up to order 3 for all points in x (an estimate of the fixed point). As a result we obtain 
a symplectic map (up to order 3) 

T : (qx,...,q n ,Px,---,Pn) -> (oi, ■ ■ ■ , On, Pi, ■ ■ ■ ,Pn)- 

We recall that from the previous step we know that the eigenvalues A = (Ai, . . . , A2 n ) of the linear part 
of T lie on the unit circle. 

As a second step we pass the linear part to a diagonal form 

A = diag(Ai, A 2 , . . . , A 2n )- 

To this end we use the linear change, B, given by a matrix formed by eigenvectors corresponding to the 
above eigenvalues. As we do not know the exact eigenvalues, the estimate of an eigenvector corresponding 
to Ai, has to be valid for all Ai £ Aj. To ensure that B is symplectic we use the following lemma and 
replace the previous eigenvectors by suitable multiples of them. 

Lemma 2. Let A G M. 2nx2n be a symplectic matrix with eigenvalues (Ai, A„, A„_|_i = Ai, . . . , X2„ = A„) 
such that \Xi\ = 1 and Ai ^ Xj for i,j = l,...,2n and i ^ j. Let (ei,...,e2 n ) be corresponding 
eigenvectors. If ef J&i+n = 1 for i = 1, . . . , n then the eigenbasis B — (ex, &2, ■ ■ ■ , s-in) is symplectic. 

Proof: To be symplectic matrix B needs to satisfy B T JB = {efJej} = J. Due to antisymmetry we 
can assume that i < j. For j = i + n we already have that ej ' Jei+ n = 1. For j ^ i + n we have AiAj ^ 1 
and therefore efJej = ej A T JAej = XiXjefjej =0. □ 
In our implementation initially £j and ei+ n are complex conjugate vectors, therefore efjei+ n = ic for 
some c £ I. We want to scale those vectors to get efJei+ n = 1 and additionally ej = iei +n . This is 
possible only if c < 0. Therefore if c > we simply interchange the indices of corresponding eigenvalues 
and eigenvectors. Finally we set 

Let u = (ux, • ■ • , U2n) be the new coordinates in this basis. The above scaling implies that Ui = i Wi+ n 
for i = 1, . . . , n. The final form of the symplectic transformation up to order 3, before starting the normal 
form computation, is 

S = B- 1 o T o B : u -> Au + S 2 + S 3 , 
where S2 and £3 denote quadratic and cubic terms of S(u) respectively. 
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Let us denote as U = {U\, . . . , U-in) the coordinates of the normal form. We achieve a normal form in 
two steps, by cancelling first all the terms of degree two in S and then the terms of degree three, except 
the unavoidable resonances. For the first step we should select, in principle, a transformation of the form 

u = N 1 (U) = U + Q{U), 

where Q(U) — (Qi, ■ • ■ , Q2n) T are quadratic terms. To cancel the terms of degree two in S we require 
that 

(4) S o JVi = ATi o A 

holds up to degree two. Let us express this in coordinates. Assume that the quadratic terms of S and 
Ni are written, respectively, as 

(5*2 )j = c i.* u *' Qi = J2 d j,kU k , for j = 1, 2n, 

|fe|=2 |fe|=2 

where k is a multiindcx. The condition Q allows to obtain 

dl ' k = x~¥ 

for all the required indices i,k. If there are no resonances of order 2, then all denominators are different 
from zero. We ensure this using rigorous estimates A^. 

However, the map N\, as it was introduced, is not a symplcctic map. Its differential satisfies the 
symplecticity conditions only to order 1 in U and we need additional cubic terms to satisfy them to order 
2. This suggests to define the symplectic transformation as the time-1 map of some Hamiltonian W 

u = N{U) = V f = r{U). 

To determine W we require 

dW dW 

ou^; = Qi> -wr Qj+n ' forj = 1 "--'" 

As the difference between the map and the identity starts with quadratic terms, the Hamiltonian starts 
with cubic terms. Then the time-1 map adds terms of degree 3 to the initial ones. Finally we obtain the 
components of N(U) as 

(N(U))\ j = U j + Q j + ±{Q j ,W}, 
where { , } denotes the Poisson bracket. This produces the normal form to order two as 

S : U ^ AU + S 3 {U), 

where 5*3 (U) are the cubic terms. 

The last step is to remove all cubic terms except those corresponding to unavoidable resonances. We 
know that there exists a symplectic, near the identity transformation that will cancel the non-resonant 
cubic terms leaving the resonant terms unchanged. Hence we can simply set to zero all terms in Ss(U) 
for which we are able to verify non-resonant condition. 

Finally we obtain the normal form to order three (we use z as a new variables) 

S:zi-> Az + T 3 (z), 

where 13(2:) are the cubic terms corresponding to unavoidable resonances. 

It remains to check the non-degeneracy condition which allows to apply KAM theorem. This is now 
easy and will be done directly on the examples. 

5. Results 

In this section we present the results of the application of the algorithm to the examples introduced 
in sections [2] and [3l The computed values are often very thin intervals with lower and upper bound 
having many common leading digits. To increase readability when printing intervals we put first those 
common digits and then the remaining digits of lower and upper bound as subscript and superscript 
correspondingly, e.g. 

123.456789^1 = [123.4567891234,123.4567895678], 
-123.456789^ = [-123.4567895678,-123.4567891234], 
-aoooooooooi = [-0.0000000001,0.0000000001]. 
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period 


3 


16 


precision 


double (52 bits) 


multiprecision (100 bits) 


X 


/ 0.35271855733^^| \ 
0.547882838S1 J 


/ +0.0000000000000001 \ 
[ -0.0000000000000001 | 

^ 0.517806657995457| J 


diam(a;) 


7-10- 12 


1Q -23 


Ai, A2 


0.5391 + i 0.841981 


-0.93239050| ± i0.36145255| 


Cl 


OQ315.580 ■ -1 o Q 345.393 
Oa 149.071 lioJ 180.278 


137393j;°H-i532621.^ 




OQ307.807 :ioq337.719 
° y 156.837 lloy 187.959 


-137393j;gi-i532621.J^ 


d 


ORQ42.269 , ; +18.283 
ZDO 05.594 ' 1 -18.281 


-234523.911 + i^ 



Table 1. Rigorous estimates of the initial condition x, eigenvalues Ai, A2, coefficients 
of normal form Ci , c 2 and torsion d for the proof of the KAM stability of orbits of period 
3 and 16 for the quartic potential. 



To be rigorous the presented intervals are rounded outwards, e.g the interval [10 -21 , 10 8 ] when rounded 
to 3 decimal places reads 0.00J. That sometimes can suggests that consecutive interval operations are 
not correct but in fact they are performed with much higher precision than the one displayed, e.g. the 
result of 1000 • 0.00^ can still be equal to 0.00 J. 

5.1. The quartic potential example. For an approximated initial condition as given in section [2] we 
carry out the algorithm described in section |4] The computed values are displayed in the Table [1] From 
the first step of the algorithm we obtain a rigorous estimate of the initial condition: x . We use them to 
get enclosures for the eigenvalues: Ai and A 2 . The Birkhoff normal form of the third and the sixteenth 
iterate of the Poincare map, respectively, computed at the fixed point x £ x in both cases is proved to be 

for some Xj £ \j and Cj € cj for j = 1,2. It is enough to work only with the first equation, because 
z\ = \Z2t Hence the map (up to order 3) reads 

(6) z\ h+ z\ exp(27ri (a + d\zi\ 2 )) 

where Ai = exp(27ria) and d = 2 C \ . Finally, for the twist map ([6]) if d (the torsion, or twist coefficient in 
that case) is different from zero then the fixed point is KAM stable (see [3j[2T]). The computed rigorous 
bounds d, shown in Table [l] verify that d is not zero for both orbits. Therefore they are KAM stable. 

The proof of the existence of the period 16 orbit using double precision did not succeed, and we were 
forced to use multiprecision interval arithmetic and Taylor method of higher order. This increases in a 
significant way the computational time. 

5.2. The figure eight orbit. We start with the coordinate system introduced in section[3]in which the 
Hamiltonian of the planar 3-body problem has the form ^fy. The variables are (<?i, <?2i <?3jPijP2jP3)- Fot 
the figure eight orbit and nearby orbits we select as Poincare section the passage through q 3 = 0, that 
is, when the three bodies are aligned. We recall that this collinear passage happens 6 times in a full 
revolution. Then, given qi,q2,pi and P2 and the value of the energy H = h we can determine > 0. 
The Poincare map defines a symplectic 4D map on the fixed level of energy and M = 0. Local variables 
which are canonically conjugate are q\,q2,p\,P2- 

Rough initial conditions in these variables if we start when the bodies are aligned, on the x axis, with 
mi to the right and 7712 at the center, are 

qi = 1.9909837697297968, q 2 = 0.9954918848648984, q 3 = 0, 
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Pi = -0.34790196497952825, p 2 = 0.69580392995905650, p 3 = 1.0678596267584018. 
Of course, p 3 can be recovered from the level of energy, which was fixed to 

h = -1.2929708570. 

The corresponding period is close to 2ir. Of course, by the homogeneity of the potential, any period or 
any negative value of h are equivalent. 

Now everything is prepared to carry out the algorithm described in section |4] The computations 
are done using multiprecision intervals with 200 bits of mantissa. For ODE integration we used Taylor 
method of order 50. As a result we obtained very sharp rigorous enclosures for the initial conditions of 
the Eight 



x = (1.990983769917805^, 0.995491884958902^, -0.34790196496310^, 0.69580392992620 



Next we computed the Taylor expansion of the Poincare map at the fixed point x G x and we proved 
that all eigenvalues of the linear part lie on the unit circle. Their rigorous estimates are: 

Ai, A 3 = 0.9985999820920383 Ti0.052896840792060g 



A 2 ,A 4 = -0.29759666751871^ ±i 0.9546916902758485 



The final normal form is 



(7) 



where 



/ Zl ^ 




( Aizi - 


h cuzfz 3 + 


C12Z1Z2Z4 \ 






A 2 z 2 - 


h C21Z1.Z2.Z3 


+ C22Z2M 






A3Z4 - 


h c 3 iziz| + 


C32Z2Z3Z4 


V * ) 




^ A4Z4 - 


1- C42Z1Z3 z 4 


+ c 4 iz 2 z| J 



en e -235.971 
on G -6.56^ + 
c 3 i 6 235.971 + 
C42 G 6.563 + 



+ il2.f 00 
121.05 
i 12 



499' 
6 

'->■ 



da G 22. 02^ 



i 1.161, 



500 
■499; 



i21.05|, 



C21 G 
C32 G - 



-0.15f 
22. 02^ 



i0.48 
i 1.16 



6) 

7 

6> 



C41 G 0.15? + i 0.48^. 



Due to the form of the eigenvalues and to the fact that no resonances appear to order 3, except the 
unavoidable ones, and using the fact that the second and fourth variables are related to the complex 
conjugates of the first and third ones, we only need to work with the two complex variables Zi,Z2- 

At this point the map reads as 

A1Z1 + cn\zi\ 2 zi + c l2 \z 2 \ 2 z 1 



(8) 







M 







Let us write Xj 



(9) 



A 2 Z 2 + C 2 l|zi|^Z 2 + C22|z 2 | 2 Z 2 

exp(i27raj) for j = 1, 2. Then the map can be written in the form 

zi exp(z27r(ai + dn|zi| 2 + rfi 2 |z 2 | 2 )) 
z 2 exp(i27r (a 2 + d 2 i \z x j 2 + d 22 \z 2 \ 2 ) ) 



7 




H 


\ Z 2 





where djk — Cjk/Xj for j = 1, 2, k = 1, 2. In the above expression we let aside, as in all the computations, 
terms of order 4 or higher. To obtain © we take first X1Z1 and A2Z2 as factors in © and compute 
logarithms. Note that the values of the coefficients djk, j = 1, 2, k = 1, 2 must be real, although from 
rigorous computations we get complex intervals around some real point. It should also hold that the 
matrix formed by the djk coefficients is symmetric, a fact which is compatible with the results obtained 
in the computations. 



For the figure eight orbit we have 

du G -37.6 



hO.001 
-0.001' 



di2 G 3.51^ + i 



-0.001 

-0.001' 



J cQtil 1 : +0.001 
"21 G d.Olg + l_ .ooi> 



d 22 G 0.08 2 + i 



+0.001 
-0.001 



For the map (|9|) the non-degeneracy KAM condition is simply that the determinant of the torsion d 
d\\d22 — ^12^21 must be different from zero. Because for the Eight we have 



+0.001 
-0.001 



ded= -15.372 s 

this finishes the proof that it is KAM stable on the level of fixed energy and angular momentum M = 0. 
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5.3. Stability of rotating Eights. Exactly the same algorithm as for the figure eight orbit can be 
used for the rotating Eights introduced in section 13.21 The obtained Birkhoff normal form and torsion 
condition are the same. Therefore, keeping the same notation as in the previous section, we only list in 
Table [2] the rigorous estimates that verify KAM stability of rotating Eights for the two selected values of 
angular momentum M x = 0.0048828125 and M 2 = 0.1484375. 



M 




0.0048828125 
/ 1.990931659347265| ^ 






0.1484375000 
f 1.951511291544930| ^ 




X 




1.008912014249292J 






1.3751748670457711 






-0.357203564790858| 
^ 0.695944468385828| ) 






-0.7065769953712011 
^ 0.835536267754800| ) 




diam(a;) 




1Q -52 






2 • 10~ 30 




An., A3 


0.99851219^ =pi 0.05452883 


5 
i 


0.65401844^ =f i0.756478ggg 


A2, A4 




0.29874539? ± i 0.95433285? 




0.620922541 ± i 0.783871921 


D = (d jk ) 


/ -36.39I + i±g:88i -3-soI + ii£88i ^ 

1. -3.501 + i±S:88i 0-08l + i±°;^ ) 


/ n c50 , : +0.001 9 ol , ■ +0.001 \ 
/ U.0 49 -I- 1_q 001 O.zo -t- i_ 001 \ 

I Q OOl 1 i +0-001 1 rrcS , : +0.001 1 
V • 3 - /0 o + 1 -0.001 1-007 + 1 -0.001 / 


d = det(D) 




1 k Aq7 1 ; +0.001 
-ia.4y 6 + i_ .ooi 






11 m 8 -i- i +o-ooi 

ll.Dl 7 + 1_ .001 





Table 2. Rigorous estimates of the initial condition x, eigenvalues Aj and torsion 
d for the proof of KAM stability of rotating Eights with angular momentum Mi — 
0.0048828125 and M 2 = 0.1484375. 
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